Scaling of Earthquake Models with Inhomogeneous Stress Dissipation 
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Natural earthquake fault systems are highly non- homogeneous. The inhomogeneities occur be- 
cause the earth is made of a variety of materials which hold and dissipate stress differently. In 
this work, we study scaling in earthquake fault models which are variations of the Olami-Feder- 
Christensen (OFC) and Rundle- Jackson-Brown (RJB) models. We use the scaling to explore 
the effect of spatial inhomogeneities due to damage and inhomogeneous stress dissipation in the 
earthquake- fault- like systems when the stress transfer range is long, but not necessarily longer than 
the length scale associated with the inhomogeneities of the system. We find that the scaling depends 
not only on the amount of damage, but also on the spatial distribution of that damage. 



PACS numbers: 

I. INTRODUCTION 

The spatial arrangement of fault inhomogeneities is de- 
pendent on the geologic history of the fault. Because 
this history is typically quite complex, the spatial dis- 
tribution of the various inhomogeneities occurs on many 
length scales. One way that the inhomogeneous nature 
of fault systems manifests itself is in the spatial patterns 
which emerge in seismicity graphs P, Q . 

Despite their inhomogeneous nature, real faults are 
often modeled as spatially homogeneous systems. One 
argument for this approach is that earthquake faults 
have long range stress transfer [3|, and if this range is 
longer than the length scales associated with the inho- 
mogeneities of the system, the dynamics of the system 
may be unaffected by the inhomogeneities. However, it 
is not clear that this is the case. Consequently it is im- 
portant to investigate the situation in which the stress 
transfer range is comparable to or less than the length 
scales associated with the damage or stress dissipation 
inhomogeneities . 

In this work, we study scaling in cellular automa- 
ton models of earthquake faults. We use a variation 
of a model introduced initially by Rundle, Jackson and 
Brown (RJB) and re-introduced independently by Olami, 
Feder and Christensen (OFC) to explore the effect of 
spatial inhomogeneities in earthquake-fault-like systems 
when stress transfer ranges are long, but not necessarily 
longer than the length scales associated with the inho- 
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mogeneities of the system 0, |5[ . For long range stress 
transfer without inhomogeneities, as well as randomly 
distributed inhomogeneities [6] such models have been 
found to produce scahng similar to Gutenberg-Richter 
scaling found in real earthquake systems [7] . It has been 
shown that the scaling found in such models is due to a 
spinodal in the limit of long range stress transfer (8|, Q . 

In the earthquake lattice models we use in this work 
we introduce inhomogeneities in the way that stress is 
dissipated. Stress is dissipated both at the lattice site of 
failure (site dissipation) and at neighboring sites which 
are damaged (damage dissipation). Spatial inhomo- 
geneities are incorporated by varying this stress dissipa- 
tion throughout the system in different spatial arrange- 
ments. We find that the scaling for damaged systems de- 
pends not only on the amount of damage, but also on the 
spatial distribution of that damage as well as the relation 
of the spatial damage or dissipation to the stress transfer 
range. Studying the effects of various spatial arrange- 
ments of site dissipation provides insights into how to 
construct a realistic model of an earthquake fault which 
is consistent with Gutenberg-Richter scaling. 



II. MODEL 

We use a two-dimensional cellular automaton model 
of an earthquake fault which is a variant of the RJB 
model 0, [ll| and closely resembles the OFC model [5[. 
We begin with a two-dimensional lattice, where each site 
is either dead (damaged) or alive (active). Each live site 
i contains an internal stress variable, cri(t), which is a 
function of time. All stress variables are initially below 
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a given threshold stress, a* and greater than or equal 
to a residual stress (both of which we assume to be 
spatially homogeneous.) Sites transfer stress to z neigh- 
bors. Neighbors are defined as all sites within the transfer 
range, R. Initially we randomly distribute stress to each 
site so that < ai < . We then increase the stress on 
all sites equally until one site reaches cr^. At this point, 
the site at the threshold stress fails. When a site fails, 
some fraction of that site's stress, given by a^(cr^ — a^^r^), 
is dissipated from the system, where ai is a parameter 
that characterizes the fraction of stress dissipated from 
site i, and is a random fiatly distributed noise. The 
stress of the site is lowered to ± and the remaining 
stress is distributed equally to the site's z neighbors. 

To model more realistic faults, we use systems which 
are damaged^ meaning they have both alive sites, which 
obey the rules out fined above, and dead sites which do 
not hold any stress. Following Serino, et al in ad- 
dition to the stress dissipation regulated by the site dis- 
sipation parameter, a^, we specify that any stress which 
is passed to a neighboring dead site also gets dissipated 
from the system. We can therefore regulate the spatial 
distribution of stress dissipation from the system with the 
distribution of the and the placement of dead sites on 
the lattice. After the initial site failure, all live neighbors 
are then checked to see if their stress has risen above a^. 
If it has, this site goes through the same failure procedure 
outfined above until all sites have stress below . The 
size of the avalanche is the number of failures that stem 
from the single initiating site. We refer to this whole 
avalanche process as a plate update. 

Because stress is dissipated from the system both at 
the site of failure (as regulated by a^) and through dead 
sites which may be placed inhomogeneously throughout 
the system, we may think of each site i as having a 
parameter which incorporates both types of dissipation, 
7i = 1 — 0i(l — Qfi), where is the fraction of live neigh- 
bors of site i. The mean value 7 = 7i/A^a, where Na is 
the number of live sites, is the average fraction of excess 
stress dissipated from the system per failed site. 

We will want to compare the scaling in these systems 
with the scaling in systems where the damage distribu- 
tion is uniform. It has been found for these OFC 
type models with no spatial inhomogeneities (homoge- 
neous damage and constant a^) that in the mean field 
limit the number of avalanche events of size s obeys the 
scaling form 



n{s) 



/s- 
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The quantity Ah, which is a function of the fraction of 
dead sites, is a measure of the distance from the spinodal 
and r = 3/2. (Note that n{s) is the number of events of 
size 5, which is the non-cumulative distribution, rather 
than the number of events of size s or smaller, which is 
the cumulative distribution often discussed in relation to 
the Gutenberg- Richter law.) We know from Ref. [12] that 
long range damaged systems with a fraction cj) of live sites 
and constant site dissipation parameter ai are equivalent 
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TABLE I: Averages and variances of ji for the distributions 
of dead sites shown in Fig. [U The total number of dead sites 
is equal to 25% of the lattice for all distributions. 



to undamaged systems with site dissipation parameter 
= 1 — — a). These systems approach the spin- 
odal {Ah 0) as the stress dissipation from the system 
vanishes: 0^1 and a ^ 0. Physically, stress dissipa- 
tion from the lattice system suppresses large avalanche 
events. 



III. QUALITATIVE BEHAVIOR OF SCALING 

First we study the case with constant ai = a and dam- 
age distributed inhomogeneously throughout the system. 
In Fig. [1] we show two dimensional lattices of linear size 
L = 256 and 25% of the sites dead. The latt ices h ave var- 
ious distributions of the dead sites. Figure l(a)| has the 
dead sites randomly distributed throughout the system. 
In the long range limit, this corresponds to homogenous 
damage studied in Ref. [12]. Figures 1(b) - 1(d) incorpo- 
rate some clustering of dead sites. Figure 1(b) has blocks 
of randomly distributed dead sites with cascading length 
scales where the linear block sizes range from 1 to L/S. 
Each block has a fraction p of randomly distributed dead 
sites, where p varies from block to block. The values 
of p are selected from a random Gaussian distribution. 
Figure p^(c)] has dead blocks with cascading length scales 
where the linear block sizes also range from 1 to L/S. 
Figure 1(d) has randomly distributed dead blocks with 
blocks of linear size of L/16 only. To characterize each 
configuration in Fig. [TJ we calculate 7, and the variance 
of 7i for an interaction range R = 16 and ai = OVi. The 
results are summarized in Table HI 

Figure [2] shows n(s), the numerical distribution of 
avalanche events of size 5, corresponding to the various 
distributions of damage in Fig. [TJ We find that the scal- 
ing behavior of systems with damage depends not only on 
the total amount of damage to the system but also on the 
spatial distribution of damage. In particular, large events 
are suppressed more for lattices with damaged sites dis- 
tributed more homogeneously. Because these lattices are 
of equal size, have the same number of damaged sites and 
the same stress transfer range(i? = 16) the differences in 
the large event behavior are not due to the finite size of 
the lattice or the finite number of active sites in the lat- 
tice. Furthermore, the results of Sec. IIII B1 indicate that 
the effect is due to the spatial distribution of 7^s and 
does not even require that the lattice be damaged. The 
calculated quantities in Table |I] would appear to indicate 
that the large event suppression is correlated with both 
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(a) random 




(c)cascading dead blocks 



(d)dead blocks 



FIG. 1: Various configurations of 25% dead sites (in black) 
for a lattice with linear size L — 256. Lattices contain (a) 
dead sites distributed randomly, (b) blocks of various sizes, 
where each block has p randomly distributed dead sites with 
p varying for each block, (c) dead blocks of various sizes, (d) 
dead blocks of a single size. 
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FIG. 2: Numerical distribution of avalanche events of size s for 
various spatial distributions of dead sites. Data corresponds 
to lattices in Fig. [T]with interaction range R = 16. 
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FIG. 3: Numerical distribution of avalanche events of size s for 
blocks of dead sites of linear size b. (Figure p, (d) | corresponds 
to R/b = 1.) The size of the system is L = 256 and the 
interaction range is = 16. The line is drawn to show that 
the data is approaching a power law with exponent —3/2. 



A. Length Scales 

For any given distribution of damage, the system will 
act as if the damage is homogeneous if the stress trans- 
fer range is long enough compared to the length scales of 
damage of the lattice. To illustrate the importance of rel- 
ative length scales, we consider the case of a single length 
scale associated with damaged areas. We place blocks of 
damaged sites of size, 6, randomly in the system which 
has constant = a. See, for example. Fig. 1(d) As we 



vary the ratio R/b^ the measured value of 7 varies from 
7 = a for R/b <C 1 to 7 = 1 - (/)(1 - a) for R/b > 1. 
In the former case, the live domains of the system ap- 
pear nearly homogeneous with (pi = 1 except near the 
boundaries of dead blocks. The latter case is the limit 
of homogeneously distributed damage. In both limiting 
cases, the variance of 7^ is small and the scaling is equiva- 
lent to the scaling for an undamaged system with a' = 7. 

In Fig. [31 we compare systems with randomly dis- 
tributed dead blocks of various length scales, R = 16, 
a = 0, and 25% total damage. As R/b gets small, the 
values of 7 also get small. The distribution of the 
corresponding data approaches a power law with the ex- 
ponent —3/2, which is the form of the distribution of a 
system at the spinodal. 



B. Spatial Distributions of Dissipation 



higher values of the average dissipation parameter 7 and 
lower values of the variance of 7. 

In order to better understand these results, we now 
study the effect of the interaction range relative to the 
length scales of inhomogeneities and the effect of cluster- 
ing of dead sites. 



The spatial distribution of damaged sites determines 
the spatial distribution of 7^ values. A more direct way to 
control the numerical and spatial distributions of 7^ is to 
use undamaged systems and vary the values of ai . In this 
way, we can isolate the effects of spatial redistribution of 

values while holding the numerical distributions of 7^ 
constant. 

We present data in Fig. [H for three systems with site 
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dissipation only; that is, they have no damage and 7^ = 
ai for each system. We see that 7 = 0.5 for ah three 
systems from the numerical distributions of ai values, 
p{ai). 

The two systems labeled "Gaussian Split" and "Gaus- 
sian Centered" both have a uniform spatial distribution 
of ai values. However, as shown in Fig. IH the values of 
ai for the "Gaussian Centered" system have a Gaussian 
distribution centered about = 0.5, while the values of 
ai for the "Gaussian Split" system have partial Gaussian 
distributions and are clustered near the values of ai = 
and ai = 1. Thus, the variance of ai values for the 
"Gaussian Centered" system is less than the variance for 
the "Gaussian Split" system. We see from the numerical 
distribution of avalanche events, n{s)^ in Fig. |4] that the 
"Gaussian Split" system has slightly larger events than 
the "Gaussian Centered" system, consistent with the ob- 
servations above that large event suppression correlates 
with low variances of a^. 

However, we see by studying the "Clustered Blocks" 
system that spatial distributions of ai have a robust ef- 
fect on scaling, even when the variances of ai are the 
same. The "Gaussian Split" and "Clustered Blocks" sys- 
tems have nearly the same numerical distributions of ai 
(Fig. HI, and therefore have the same value of the vari- 
ance of ai. The spatial distributions of these two cases, 
however, are different: the "Gaussian Split" system has a 
uniform spatial distribution of ai values, while the "Clus- 
tered Blocks" system has high (and low) ai values clus- 
tered together into blocks as shown in the inset of Fig. |4j 
Despite having equal values of Oi and equal variances 
of ai values, the "Clustered Blocks" system experiences 
much larger events (by an order of magnitude). 

Evidently, the larger events depend crucially on the 
spatial clustering of low dissipation sites. This is because 
failing sites with low values of 7^ pass along a high per- 
centage of excess stress, encouraging the failure of neigh- 
boring sites. Thus, a large earthquake event is more likely 
to occur if the initial site of failure is well connected to 
a large number of sites with low dissipation parameters. 
In our system, connectedness is determined by spatial lo- 
cality, so we require large clumps of sites with low values 
of 7i in order to allow for the occasional large earthquake 
event. 



IV. GUTENBERG-RICHTER SCALING 

The Gutenberg-Richter scahng law states that the cu- 
mulative distribution of earthquake sizes is exponential in 
the magnitude [7]. In terms of the seismic moment, which 
has succeeded the Richter magnitude as the appropriate 
measure for earthquake sizes, the law may be reframed 
to state that the cumulative distribution of earthquake 
sizes, Nm^ is a power law in the seismic moment, M [6]. 

Nm - M-l" , with 13= J, (2) 
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FIG. 4: Comparison of three lattice systems with no damage 
and two different distributions of p{ai) shown in the top fig- 
ure. The lattices labeled by "Gaussian Centered" and "Gaus- 
sian Split" are distributed uniformly in space, while the spa- 
tial distribution of "Clustered Blocks" is shown in the in- 
set. The bottom plot shows the numerical distribution of 
avalanche events of size s. 

and b is the so-called b value of the Gutenberg-Richter 
law which has been measured for many real earthquake 
systems. The seismic moment M is proportional to the 
size of the earthquake in this model jl3|- Therefore, the 
relation appropriate for the systems considered in this 
work is the cumulative distribution of earthquake size: 

Ns ^ s-^, (3) 

or the corresponding non-cumulative distribution 

Us ^ , with f = /3 + 1. (4) 

Serino et al [6] construct a model for an earthquake 
fault system consisting of an aggregate of lattice models, 
where each lattice has a fraction q of homogeneously dis- 
tributed dead sites and q varies from to 1. The weight- 
ing factor Dq gives the fraction of lattices with damage 
q. Considering the weighting factor to be constant with 
all values of q contributing equally to the fault system, 
they find a value of f = 2. They also consider a power 
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law distribution of Dq and fit the exponent to correspond 
to Gutenburg-Richter b values found in real earthquake 
systems. 

There are two important differences between the model 
considered by Serino et al and our work: In the model 
treated by Serino et al 

1. The damage is distributed homogeneously. 

2. The individual lattices with homogeneous damage 
q are non-interacting. 

We investigate both the effect of the spatial arrangement 
of the damage and its relation to the stress transfer range 
as well as the effect of stress transfer between regions with 
different levels of damage. 

We construct two lattice systems with a uniform nu- 
merical distribution of 7^ which have scaling consistent 
with Serino et. al.'s systems with constant distribution 
Dq. The first model essentially pieces together many ho- 
mogeneous lattice systems: The numerical distribution 
of ai values is uniform between and 1, but spat ially 
arranged into Nb blocks of linear size B (see Fig. |5(a 
inset), where each block contains a random distribution 
of ai values within an interval of size 1/Nb- There are no 
dead sites, so that = 7^. The effects of the boundaries 
between the blocks should be negligible \i B <^ R. In 
Fig. 5(a)[ we present data from a system with L = 512, 

= 16, and B = 64. The straight line shows the best fit 
to a power law with exponent f o:^ 2.07, which is consis- 
tent with the results for the aggregate lattice system of 
Serino et al with Dq = 1. 

We find that the size of the blocks, 5, need not be 
the same for different values of a^. It is important that 
the boxes with lower values of 7 be large enough to ac- 
commodate large avalanche events, but blocks with large 
ai may be small because they are more likely to seed 
small avalanches. With this in mind we construct a lat- 
tice system with cascading length scales of blocks where 
the largest blocks have the lowest ai values and decreas- 
ing sized blocks have increasing values of ai. The scaling 
results are shown in Fig. |5(b)[ with a best fit power law 
with exponent f o:^ 2.04. 



V. CONCLUSIONS 

We have studied both damage and site dissipation to 
inform the development of models of realistic earthquake 
faults with inhomogeneous stress dissipation. Spatially 
rearranging dead sites on a given lattice affects the nu- 
merical distributions of the effective stress dissipation 
parameters and the scaling behavior of large avalanche 
events, depending on the homogeneity of the damage 
and the length scales associated with the clustered dead 
sites. However, by studying site dissipation we find that 
spatial distributions of dissipation parameters crucially 



affect scaling behavior even when the numerical distri- 
butions of dissipation parameters are the same. 
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FIG. 5: Numerical distribution for avalanche events of size s 
for systems with uniform numerical distributions for ai, but 
non-uniform spatial distributions of ai which are shown in 
the insets. Slopes of best fit lines in red are a) f ^ 2.07 and 
b) f ^ 2.04. 



Sites with lower stress dissipation, even if only partially 
distributed throughout the lattice but clumped together, 
allow for larger avalanche events. We have found models 
for earthquake fault systems which have avalanche event 
size scaling which is consistent with the new paradigm 
for Gutenberg-Richter scaling proposed by Serino et al. 
The models studied here go beyond those previously pro- 
posed by incorporating inhomogeneities into the lattice 
and allowing areas with different characteristic dissipa- 
tion rates to interact. 
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